% Band structure of InSb following Kane's 8x8 k.p Hamiltonian including non-parabolicity 
% corrections, calculated with different approximations. 
%
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM.
% Date: 18 October 2023
% Required package: None
%
% References:
% E. O. Kane, Band structure of indium antimonide, J. Phys. Chem. Solids 1, 249-261 (1957).
% E. O. Kane, The k.p method, in Semiconductors and Semimetals, edited by Willardson and Beer, vol. 1, 75-100 (1966). 
% Stefan Zollner and Jose Menendez, Nonparabolicity and direct gap absorption of InSb
% (unpublished, 2023). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Physical constants (energy units in eV unless noted) from CODATA at NIST
hbar=6.582119569E-16; % eV*s
e=1.602176634E-19; % A*s
m0=9.1093837015E-31; % kg
eps0=8.8541878128E-12; % F/m
kB=0.08617333262E-3; % eV/K Boltzmann constant
c=299792458; % m/s speed of light in vacuum
%
kBJoule=kB*e;
h=hbar*2*pi; % eV*s
hJoule=h*e; % J*s
hbarJoule=hbar*e;
hbarc=hbar*c*100; % eV*cm
hbar2m0=hbar*hbarJoule/m0*1E20; % hbar^2/m0 in Angstrom^2
Rydberg=e^4*m0/(8*eps0^2*hJoule^2)/e; % eV, hydrogen ground state energy
Hartree=2*Rydberg; % eV, atomic unit for energy 
aBohr=(hJoule^2*eps0)/(pi*e^2*m0)*1E10; % Angstrom
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Material parameters for InSb (compare Ioffe Institute web site)
a=6.47937; % lattice constant of InSb in Angstrom at 298 K, Landolt-Boernstein vol. 17a
E0=0.237; % eV, from Lawaetz 1971, at 0 K
Delta0=0.81; % eV, spin-orbit splitting
EP=23.1; % eV, momentum matrix energy, from Lawaetz 1971, at 0 K
% EP=0.44*Hartree*2; % convert 2P^2 from atomic units to eV, from Kane 1957
% effective mass from perturbative solution of cubic equation
me=1/(1+EP/3*(2/E0+1/(E0+Delta0))) % should be about 0.0136 (Lawaetz)
% effective mass in the limit of large SO splitting
% me=3*E0/(2*EP+3*E0)
mhh=0.47; % from Lawaetz 1971
% effective light hole mass in the limit of large SO splitting (same as perturbative)
mlh=3*E0/(2*EP-3*E0)
% mlh=0.015; % from Lawaetz 1971
mso=0.107; % from Lawaetz 1971
% Calculated parameters
muhh=mhh*me/(mhh+me); % reduced (optical) mass, compare Emminger 2022, Table II
mulh=mlh*me/(mlh+me)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate dispersions of valence and conduction bands for InSb.
% Use approximation for large spin-orbit splittings with Taylor series expansion.
Ehh=@(k)hbar2m0*k.^2/2;
Ecp=@(k)E0+hbar2m0*k.^2/(2*me);
Ec=@(k)hbar2m0*k.^2/2+E0/2*(1+sqrt(1+8/3*hbar2m0*k.^2/2*EP/E0^2));
Ec1=@(k)E0+hbar2m0*k.^2/2 .* ...
    (1+2/3*EP/E0-4/9*hbar2m0*k.^2/2*EP^2/E0^3);
Ec2=@(k)E0+hbar2m0*k.^2/2 .* ...
    (1+2/3*EP/E0-4/9*hbar2m0*k.^2/2*EP^2/E0^3+2*8/27*hbar2m0^2/4*k.^4*EP^3/E0^5);
%
Elhp=@(k)-hbar2m0*k.^2/(2*mlh);
Elh=@(k)hbar2m0*k.^2/2+E0/2*(1-sqrt(1+8/3*hbar2m0*k.^2/2*EP/E0^2));
Elh1=@(k)hbar2m0*k.^2/2 .* ...
    (1-2/3*EP/E0+4/9*hbar2m0*k.^2/2*EP^2/E0^3);
Elh2=@(k)hbar2m0*k.^2/2 .* ...
    (1-2/3*EP/E0+4/9*hbar2m0*k.^2/2*EP^2/E0^3-2*8/27*hbar2m0^2/4*k.^4*EP^3/E0^5);    
%
k=-0.05:0.001:0.05;
kau=k*aBohr; % wave vector k in atomic units
% Draw a figure of energies versus k in atomic units.
figure('Name','Nonparabolicity of bands in InSb from 4x4 k.p model with large SO', ...
    'Color','white');
pl=plot(kau,Ec(k),'r-',kau,Ecp(k),'r:',kau,Ec1(k),'r--',kau,Ec2(k),'r-.', ...
        kau,Ehh(k),'b-', ...
        kau,Elh(k),'g-',kau,Elhp(k),'g:',kau,Elh1(k),'g--',kau,Elh2(k),'g-.'); 
xlabel('wave vector k (atomic units)'); ylabel('Energy (eV)');
ylim([-0.2 0.5]); xlim([-0.015 0.015]); 
pl(1).LineWidth=3; pl(5).LineWidth=3; pl(6).LineWidth=3;
legend('"exact" large SO approximation','parabolic','quartic','sextic', ...
    'Location','north'); 
legend('boxoff');
% Same figure with k^2 along the horizontal axis.
k=0:0.001:0.2;
kau=k*aBohr; % wave vector k in atomic units
% Compare Kane 1957, Figure 1
figure('Name','Nonparabolicity of bands in InSb from 4x4 k.p model with large SO', ...
    'Color','white');
pl=plot(kau.^2,Ec(k),'r-',kau.^2,Ecp(k),'r:',kau.^2,Ec1(k),'r--',kau.^2,Ec2(k),'r-.', ...
        kau.^2,Ehh(k),'b-', ...
        kau.^2,Elh(k),'g-',kau.^2,Elhp(k),'g:',kau.^2,Elh1(k),'g--',kau.^2,Elh2(k),'g-.'); 
xlabel('wave vector k^2 (atomic units)'); ylabel('Energy (eV)');
xlim([0,0.002]); ylim([-0.6,1]);
pl(1).LineWidth=3; pl(5).LineWidth=3; pl(6).LineWidth=3;
legend('"exact" large SO approximation','parabolic','quartic','sextic', ...
    'Location','east'); 
legend('boxoff');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eenergy dispersions from Vieta's solution to the cubic characteristic equation
k=0:0.001:0.1; % wave vector k for band structure plot
kau=k*aBohr; % wave vector k in atomic units
%
Asqu=@(k)hbar2m0*k.^2/2*EP/3;
a=1;
b=Delta0-E0;
c=@(k)-E0*Delta0-3*Asqu(k);
d=@(k)-2*Asqu(k)*Delta0;
% use a=1
p=@(k)c(k)-b^2/3;
q=@(k)2*b^3/27-b*c(k)/3+d(k);
%
t0=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 );
t1=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 -2*pi/3);
t2=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 -4*pi/3);
% 
Ene0=@(k)t0(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
Ene1=@(k)t1(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
Ene2=@(k)t2(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
%
figure('Name','Bands in InSb from 4x4 k.p model with Vietas equation', ...
    'Color','white');
pl=plot(kau,Ene0(k),'r',kau,Ene1(k),'g',kau,Ene2(k),'k', ...
    kau,Ec(k),'r:',kau,Elh(k),'g:',kau,Ehh(k),'b');
xlabel('wave vector k (atomic units)'); ylabel('Energy (eV)');
xlim([0,0.05]);
legend('exact 8x8 k.p Hamiltonian','','','large SO approximation', ...
    'Location','best'); 
legend('boxoff');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonparabolicity parameters calculated in the limit of large SO splittings
% alphae=(2/3*EP)^2/E0/(2/3*EP+E0)^2;
% betae=-2*(2/3*EP)^3/E0/(2/3*EP+E0)^4;
% alphalh=(2/3*EP)^2/E0/(2/3*EP-E0)^2;
% betalh=-2*(2/3*EP)^3/E0/(2/3*EP-E0)^4;
% Nonparabolicity parameters calculated from reduced masses for large SO splitting
alphaep=(1-me)^2;
betaep=-2*me*(1-me)^3;
alphalhp=(1+mlh)^2;
betalhp=2*mlh*(1+mlh)^3;
% Nonparabolicity parameters from Bartoli (1983)
% delta=Delta0/E0
% mlh=me*(delta+1.5)/(delta+1-me*(2*delta+2.5)); % same as above for perturbative
% alphaep=alphaep*(1.5+2*delta+delta^2)/(1.5+delta)/(1+delta)
% betaep=-2*(1-me)^3/(1.5+delta)^2/(1+delta)* ...
%    ( 0.25*delta^2 + me*(1.5+2*delta+delta^2)^2/(1+delta) )
% alphalhp=alphalhp*(1+0.5/delta)
% betalhp=(1+mlh)^3* ...
%    ( 0.5/delta + 0.75/delta^2 + 2*mlh*(1+0.5/delta)^2 )
%
alphae=alphaep/E0;
betae=betaep/E0^2;
alphalh=alphalhp/E0;
betalh=betalhp/E0^2;
%
meDOS1=@(eps)me*nthroot(1+alphae*eps,3).*(1+2*alphae*eps).^(2/3);
meDOS=@(eps)me*nthroot(1+alphae*eps+betae*eps.^2,3).*(1+2*alphae*eps+3*betae*eps.^2).^(2/3);
mlhDOS1=@(eps)mlh*nthroot(1+alphalh*eps,3).*(1+2*alphalh*eps).^(2/3);
mlhDOS=@(eps)mlh*nthroot(1+alphalh*eps+betalh*eps.^2,3).*(1+2*alphalh*eps+3*betalh*eps.^2).^(2/3);
%
eps=0:0.001:0.5;
figure('Name','Conduction band nonparabolicity in InSb', ...
    'Color','white');
plot(eps,meDOS1(eps),'k:',eps,meDOS(eps),'k-',eps,mlhDOS1(eps),'r:',eps,mlhDOS(eps),'r-');
ylabel('electron or hole DOS mass (m_0)'); 
xlabel('Energy below VB maximum or  above CB minimum (eV)');
legend('electron quadratic','electron cubic', ...
    'light hole quadratic','light hole cubic', ...
    'Location','northwest'); 
legend('boxoff');
% xlim([0 0.1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw a figure of energies versus k in atomic units
% with the perturbative solution to Vieta's equation. 
k=-0.05:0.001:0.05;
kau=k*aBohr; % wave vector k in atomic units
Elhpert=@(k)hbar2m0*k.^2/2 .* ...
    (1-2/3*EP/E0+ ...
    2/9*hbar2m0*k.^2/2*EP^2/(E0^2*Delta0)*(1+2*Delta0/E0));
Ecpert=@(k)E0+hbar2m0*k.^2/2/me ...
    -(1/me-1)^2*(1/(E0+Delta0)+2*Delta0/E0/(3*E0+2*Delta0))*hbar2m0^2/4.*k.^4;
figure('Name','Nonparabolicity of bands in InSb from 4x4 k.p model with SO', ...
    'Color','white');
pl=plot(kau,Ec(k),'r-',kau,Ecp(k),'r:',kau,Ec1(k),'r--',kau,Ec2(k),'r-.', ...
        kau,Ehh(k),'b-', ...
        kau,Elh(k),'g-',kau,Elhp(k),'g:',kau,Elh1(k),'g--',kau,Elh2(k),'g-.', ...
        kau,Elhpert(k),'g--',kau,Ecpert(k),'r--'); 
xlabel('wave vector k (atomic units)'); ylabel('Energy (eV)');
ylim([-0.2 0.5]); xlim([-0.015 0.015]); 
pl(1).LineWidth=3; pl(5).LineWidth=3; pl(6).LineWidth=3;
pl(10).LineWidth=1.5; pl(11).LineWidth=1.5;
legend('"exact" large SO approximation','parabolic','quartic','sextic', ...
    '','','','','','','perturbative (quartic)', ...
    'Location','north'); 
legend('boxoff');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

